function updateLTS(timeStep)
%updateLTS solves the diffEq for t = timeStep+dt and returns it
%   Detailed explanation goes here

global t;
%global SBCells;
%global RSCells;
global LTSCells;
global IB_ACells;
global Synapses;
global LTSIndices;
global IBIndices;
global Gates;
global J_SF;

dt = t(timeStep+1) - t(timeStep);


%V_aL = IB_ACells(timeStep,1,:);
%V_aL = reshape(V_aL,[size(IB_ACells,3) 1]);

for i = 1:size(LTSCells,3)
    in_s = Synapses(timeStep, :, LTSIndices(i));
    in_gates = Gates(:, LTSIndices(i))';
    
    [k1, s1] = feval(@LTS, timeStep, t(timeStep), LTSCells(timeStep,:,i),in_s, in_gates, J_SF);
    % note - caromk 5/24/2010 - these k and s updates are to normalize the
    % RK4 method with the cleaned up code, this code should not continue
    % to be built upon
    k1 = dt*k1; s1 = dt*s1;
    [k2, s2] = feval(@LTS, timeStep, t(timeStep) + dt/2, LTSCells(timeStep,:,i) + (k1)'/2,in_s + s1/2, in_gates, J_SF);
    k2 = dt*k2; s2 = dt*s2;
    [k3, s3] = feval(@LTS, timeStep, t(timeStep) + dt/2, LTSCells(timeStep,:,i)+ (k2)'/2,in_s + s2/2, in_gates, J_SF);
    k3 = dt*k3; s3 = dt*s3;
    [k4, s4] = feval(@LTS, timeStep, t(timeStep) + dt, LTSCells(timeStep,:,i) + k3',in_s + s3, in_gates, J_SF);
    k4 = dt*k4; s4 = dt*s4;
    LTSCells(timeStep+1,:,i) = LTSCells(timeStep,:,i) + (k1' + 2*k2' + 2*k3' + k4')/6;
    
    %    if t(timeStep)>50 && t(timeStep)<51
    %        LTSCells(timeStep+1,1,i) = 45;
    %    end
    
    Synapses(timeStep+1, :, LTSIndices(i)) = in_s + (s1 + 2*s2 + 2*s3 + s4)/6;
    
end

end